#install.packages("pacman")#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)#devtools::install_github("daniel1noble/orchaRd", force = TRUE)library(tidyverse)library(here)library(lme4)library(orchaRd)#library(gptstudio)library(metafor)library(patchwork)library(alluvial)library(ggalluvial)library(easyalluvial)library(ape)library(clubSandwich)library(emmeans)library(MuMIn)library(kableExtra)library(pander)# making metafor talk to MuMIneval(metafor:::.MuMIn)# install.packages("pak")#pak::pak("MichelNivard/gptstudio")
2 Getting data loaded
Code
#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))#dat_full <- read.csv(here("data/dat_28_06_2023.csv"))dat_full <-read.csv(here("data/dat_19_07_2023.csv"))# add phylogenetic tree - only topologies# TODO? - we could get better tree from birdtree.org# we can do 50 different trees as in # https://academic.oup.com/sysbio/article/68/4/632/5267840tree_top <-read.tree(here("R/birds_MA.tre"))# tree with branch lengthstree <-compute.brlen(tree_top)#plot(tree)# turning it into a correlation matrixcor_tree <-vcv(tree,corr=T)
3 Custom functions
Code
# custom functions#' Title: Contrast name generator#'#' @param name: a vector of character stringscont_gen <-function(name) { combination <-combn(name, 2) name_dat <-t(combination) names <-paste(name_dat[, 1], name_dat[, 2], sep ="-")return(names)}#' @title get_pred1: intercept-less model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred1 <-function (model, mod =" ") { name <-firstup(as.character(stringr::str_replace(row.names(model$beta), mod, ""))) len <-length(name)if (len !=1) { newdata <-matrix(NA, ncol = len, nrow = len)for (i in1:len) { pos <-which(model$X[, i] ==1)[[1]] newdata[, i] <- model$X[pos, ] } pred <- metafor::predict.rma(model, newmods = newdata) }else { pred <- metafor::predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <-tibble(name =factor(name, levels = name, labels = name), estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = model$pval,lowerPR = lowerPR, upperPR = upperPR)}#' @title get_pred2: normal model#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)#' @param model: rma.mv object #' @param mod: the name of a moderator get_pred2 <-function (model, mod =" ") { name <-as.factor(str_replace(row.names(model$beta), paste0("relevel", "\\(", mod,", ref = name", "\\)"),"")) len <-length(name)if(len !=1){ newdata <-diag(len) pred <-predict.rma(model, intercept =FALSE, newmods = newdata[,-1]) }else { pred <-predict.rma(model) } estimate <- pred$pred lowerCL <- pred$ci.lb upperCL <- pred$ci.ub lowerPR <- pred$cr.lb upperPR <- pred$cr.ub table <-tibble(name =factor(name, levels = name, labels = name), estimate = estimate,lowerCL = lowerCL, upperCL = upperCL,pval = model$pval,lowerPR = lowerPR, upperPR = upperPR)}#' @title mr_results#' @description Function to put results of meta-regression and its contrasts#' @param res1: data frame 1#' @param res1: data frame 2mr_results <-function(res1, res2) { restuls <-tibble(`Fixed effect`=c(as.character(res1$name), cont_gen(res1$name)),Estimate =c(res1$estimate, res2$estimate),`Lower CI [0.025]`=c(res1$lowerCL, res2$lowerCL),`Upper CI [0.975]`=c(res1$upperCL, res2$upperCL),`P value`=c(res1$pval, res2$pval),`Lower PI [0.025]`=c(res1$lowerPR, res2$lowerPR),`Upper PI [0.975]`=c(res1$upperPR, res2$upperPR), )}#' @title all_models#' @description Function to take all possible models and get their results#' @param model: intercept-less model#' @param mod: the name of a moderator all_models <-function(model, mod =" ", type ="homo") {# getting the level names out level_names <-levels(factor(model$data[[mod]])) dat2 <- model$data VCV1 <-vcalc(vi = dat2$Vd,cluster = dat2$SubjectID,rho =0.5)#model$data[[mod]] <- factor(model$data[[mod]], ordered = FALSE)# meta-regression: contrasts # helper function to run metafor meta-regression run_rma1 <-function(name) {rma.mv(yi = SMD, V = VCV1, mods =~relevel(dat2[[mod]], ref = name), random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),test ="t",method ="REML",sparse =TRUE,data = dat2) } run_rma2 <-function(name) {rma.mv(yi = SMD, V = VCV1, mods =~relevel(dat2[[mod]], ref = name), random =list(~1| FocalSpL,~1| RecNo,~gsub("\"", "", mod) | Obs_ID),test ="t",method ="REML",sparse =TRUE,data = dat2) }# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])if (type =="homo"){ model_all <-map(level_names[-length(level_names)], run_rma1) } else { model_all <-map(level_names[-length(level_names)], run_rma2) }# getting estimates from intercept-less models (means for all the groups) res1 <-get_pred1(model, mod = mod)# getting estiamtes from all contrast models res2_pre <-map(model_all, ~get_pred2(.x, mod = mod))# a list of the numbers to take out unnecessary contrasts contra_list <-Map(seq, from=1, to=1:(length(level_names) -1)) res2 <-map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) # creating a table res_tab <-mr_results(res1, res2) %>%kable("html", digits =3) %>%kable_styling("striped", position ="left") %>%scroll_box(width ="100%")# results res_tab}# test #all_models(mod1, mod = "Treat_mod")#all_models(mod1s, mod = "Treat_mod")
4 Preparing data set
Code
# function for calculating varianceVd_func <-function(d, n1, n2, design, r =0.5){# independent designif(design =="among"){ var <- (n1 + n2) / (n1*n2) + d^2/ (2* (n1 + n2 -2)) # variance } else { # dependent design var <-2*(1-r) / n1 + d^2/ (2*(n1 -1)) # variance } var # return variance}# getting Hedges' g - get small size bias corrected effect sizedat_full$SMD <- dat_full$d / (1-3/(4* (dat_full$NTreat + dat_full$Ncontrol) -9))# flipping d dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection# calucating Vddat_full$Vd <-with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))# extra useful function# function for getting mean and sd from median, quartiles and sample size# get_mean_sd <- function(median, q1, q3, n){# sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd# mean <- (median + q1 + q3)/3 # mean# c(mean, sd)# }# observation iddat_full$Obs_ID <-1:nrow(dat_full)dat_full$Phylo <-gsub(" ", "_", dat_full$FocalSpL)# filtering very large variance and also very small sample sizedat_int <- dat_full %>%filter(Vd <10& Ncontrol >2& NTreat >2)#dim(dat_full)#dim(dat_int)# sorting out modality stuff# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) dat_int %>%mutate(Treat_mod =case_when(Treatment =="A"~"A", Treatment =="AV"~"AV", Treatment =="AVG"~"AV", Treatment =="AVM"~"AV", Treatment =="L"~"AVO", Treatment =="O"~"O", Treatment =="OV"~"OV", Treatment =="V"~"V", Treatment =="VG"~"V", Treatment =="VM"~"V", Treatment =="VP"~"V"),# into how manyTreat_No =case_when(Treatment =="A"~1, Treatment =="AV"~2, Treatment =="AVG"~2, Treatment =="AVM"~2, Treatment =="L"~3, Treatment =="O"~1, Treatment =="OV"~2, Treatment =="V"~1, Treatment =="VG"~1, Treatment =="VM"~1, Treatment =="VP"~1),# des it have some add-onsAdd_on =case_when(Treatment =="A"~"No", Treatment =="AV"~"No", Treatment =="AVG"~"Yes", Treatment =="AVM"~"Yes", Treatment =="L"~"No", Treatment =="O"~"No", Treatment =="OV"~"No", Treatment =="V"~"No", Treatment =="VG"~"Yes", Treatment =="VM"~"Yes", Treatment =="VP"~"Yes"), ) -> dat# creating data just for A, V, and AV dat_short <- dat %>%filter(Treat_mod =="A"| Treat_mod =="V"| Treat_mod =="AV")# for add-on, we only need V and AVdat_short_add <- dat %>%filter(Treat_mod =="AV"| Treat_mod =="V")dat <- dat %>%mutate_if(is.character, as.factor)# reordering factors for better visualisation# dat$Treat_mod <- factor(dat$Treat_mod, levels = c("A", "V", "AV", "O", "OV", "AVO"))
5 Exploratory visualization
For Treat_mod (Treatment), we will only visualise A, V, and AV as O $r $, OV, and AVO are much rarer. But for Type (Trait type), we will use all data.
orchard_plot(mod1s2, mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
7.3 Comparing models
Code
# comparision modelsanova(mod1s, mod1s2)
df AIC BIC AICc logLik LRT pval QE
Full 8 1692.5188 1727.6541 1692.7637 -838.2594 NA
Reduced 6 1707.4108 1733.7623 1707.5532 -847.7054 18.8921 <.0001 NA
Code
# the effect of additions# this is a part of sensitivity analysisVCVs2 <-vcalc(vi = dat_short_add$Vd,cluster = dat_short_add$SubjectID,rho =0.5)mod5 <-rma.mv(yi = SMD, V = VCVs2, random =list(~1|FocalSpL , ~1| RecNo, ~ Treat_mod | Obs_ID), mod =~ Treat_mod*Add_on -1,test ="t",struct ="DIAG",method ="REML", sparse =TRUE,data = dat_short_add)summary(mod5)
# testing the number of stimulimod4 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Treat_No, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod4)
bubble_plot(mod4,mod ="Treat_No",group ="RecNo",xlab ="The number of simuli",g =TRUE)
Code
# Type of responsesmod2 <-rma.mv(yi = SMD, V = VCV, random =list(~1| FocalSpL,~1| RecNo,~1| Obs_ID),mod =~ Type -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod2)
# heteroscadasticity model better than the homoscedasticity model# note LifeHistory has lowest variation but this may be expected? # as it is less flexiable (e.g. the number of eggs?)anova(mod2, mod2b)
df AIC BIC AICc logLik LRT pval QE
Full 8 1757.5731 1793.2272 1757.8024 -870.7865 NA
Reduced 6 1785.4590 1812.1996 1785.5923 -886.7295 31.8859 <.0001 NA
Code
orchard_plot(mod2b, mod ="Type",group ="RecNo",xlab ="Standardised mean differnece (SMD)",branch.size =3)
7.4 Trait categories
Code
# Category of responsesmod3 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo, ~1| Obs_ID), mod =~ Category -1, test ="t",method ="REML", sparse =TRUE,data = dat)summary(mod3)
orchard_plot(mod3, mod ="Category",group ="RecNo", xlab ="Standardised mean differnece (SMD)",angle =45,branch.size =3)
7.5 Predactor guild
Code
# Predactor guild# quite heterogeneous# TODO this could be in random effects - think abou thtis a bit latermod6 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ PredGuild -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod6)
orchard_plot(mod11,mod ="ControlType",group ="RecNo",xlab ="Standardised mean differnece (SMD)")
Code
# sex# TODO - this could be interesting# what is in males and femalesmod12 <-rma.mv(yi = SMD, V = VCV, random =list(~1|FocalSpL , ~1| RecNo,~1| Obs_ID),mods =~ Sex -1,test ="t",method ="REML",sparse =TRUE,data = dat)summary(mod12)
bubble_plot(mod15,mod ="ln_duration",group ="RecNo",xlab ="log(duration in days)",g =TRUE) +geom_point(data = dat,aes(x = ln_duration, y = SMD,color = Type,fill = Type,size =1/sqrt(Vd)), alpha =0.6) +scale_color_discrete() +#+ # how to put the legend for colourguides(color ="legend")
orchard_plot(mod_full,mod ="Type",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
orchard_plot(mod_full,mod ="Treat_mod",group ="RecNo", xlab ="Standardised mean differnece (SMD)",branch.size =3)
Code
int_type <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Type")bubble_plot(int_type, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)
Code
int_trt <-mod_results(mod_full, mod ="sln_duration", group ="RecNo", weights ="prop",by ="Treat_mod")bubble_plot(int_trt, group ="RecNo", mod ="sln_duration", xlab ="ln(duration in days)",legend.pos ="top.left", condition.nrow =3)
# displays delta AICc <2candidates_aic2 <-subset(candidates, delta <5) # model averagingmr_averaged_aic2 <-summary(model.avg(candidates, delta <5)) # relative importance of each predictor for all the modelsimportance <-sw(candidates)
dat_fulle <-qdrg(object = mod_fulle, data = dat_short)# marginalized overall mean at vi = 0 and year.c = 0; also weights = "prop" or "cells" average things over proportionally. if not specified, all groups (levels) get the same weights# res_fulle1 <- emmeans(dat_fulle, # specs = ~ sqeffectN,# df = mod_fulle$ddf, # weights = "prop")
10 R Session Informtion
Code
# pander for making it look nicersessionInfo() %>%pander()